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Summary 

This  grant  supported  three  years  of  research  by  Prof.  Trefethen  and  several  graduate 
students,  including  the  thesis  research  of  Louis  Howell  (PhD  January  1990,  M.I.T.). 
Progress  was  made  in  the  following  areas: 

Algorithms  for  numerical  conformal  mapping 

•  Development  of  an  algorithm  for  computing  ideal  2D  jet  flows 

•  Development  of  an  algorithm  for  mapping  highly  elongated  polygons 

•  Comparison  of  quadrature  methods  in  numerical  conformal  mapping 

•  Proof  of  divergence  of  existing  Schwarz- Christoffel  algorithms 

•  Improvement  of  algorithms  for  mapping  circular  polygons 

•  Derivation  of  a  formula  for  mapping  highly  elongated  circular  polygons 

•  Review  of  “Schwarz-Christoffel  Mapping  in  the  1980’s” 

•  Release  of  second  edition  of  SCPACK  User’s  Guide 

Applications  of  complex  analysis  in  scientific  computing 

•  Comparison  of  matrix  iterative  algorithms  based  on  complex  approximation 

•  Discovery  that  eigenvalue-based  algorithms  are  unreliable 

•  Analysis  of  spectra  and  pseudo-spectra  of  Toeplitz  matrices 

•  Development  of  a  new  hybrid  algorithm  for  nonsymmetric  matrix  iterations 

•  Development  of  a  Lax-stability  criterion  based  on  complex  stability  regions 

•  Development  of  algorithms  for  polynomial  interpolation  and  least-squares 
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1.  Introduction 

Conformal  mapping  is  as  old  as  complex  analysis,  and  has  always  been  one  of  the 
branches  of  complex  analysis  most  closely  tied  to  applications.  Mathematically,  the 
problem  is  to  find  an  analytic  function  that  maps  a  complicated  domain  in  the  complex 
plane  onto  a  simple  one.  The  applications  come  about  because  such  maps  preserve 
solutions  to  Laplace’s  equation,  as  well  as  having  other  useful  properties,  with  the  rccv.lt 
mat  they  simplify  many  problems  in  heat  conduction,  electrostatics,  electromagnetics, 
fluid  mechanics,  probability,  and  other  fields. 

The  difficulty  with  conformal  mapping  is  that  except  in  certain  special  cases,  the  func¬ 
tion  in  question  cannot  be  determined  analytically.  One  must  resort  to  numerical  con¬ 
formal  mapping,  and  in  the  past  few  decades  a  branch  of  numerical  analysis  has  sprung 
up  with  this  name.  Numerical  conformal  mapping  is  not  as  central  a  topic  in  scientific 
computing  as  solution  of  differential  equations,  zerofinding,  or  linear  programming,  to 
name  a  few,  but  it  persists  year  after  year  as  a  sideline  of  genuine  importance. 

I  have  been  working  on  numerical  conformal  mapping  since  1978,  when  Peter  Henrici 
visited  Stanford  while  I  was  in  my  second  year  as  a  graduate  student.  My  interests 
have  centered  on  ideas  related  to  the  Schwarz-Christoffel  transformation,  which  is  a 
formula  for  the  conformal  mapping  of  polygons.  Like  all  conformal  maps,  Schwarz- 
Christoffel  maps  require  numerical  computation  in  all  but  the  most  trivial  cases.  In 
1983  I  produced  a  Fortran  package  called  SCPACK  for  carrying  out  these  computations, 
and  in  the  years  since  then,  SCPACK  has  become  widely  used. 

The  purpose  of  this  grant  was  to  support  continued  research  in  numerical  conformal 
mapping,  applications  of  conformal  maps,  and  related  topics.  Understandably  enough, 
the  starting  point  of  my  work  turned  out  to  v  ~  tV>e  Schwarz-Christoffel  transformation, 
and  my  student  Louis  Howell  and  I  made  s  ibsf<mtial  progress  on  a  number  of  prob¬ 
lems  in  this  area.  This  work  is  summarized  in  _  don  5.  Then,  while  progress  in  those 
topics  continued,  the  research  also  took  another  turn.  One  application  of  conformal 
mapping,  investigated  by  Varga  and  Reichel  and  Tal-Ezer  and  others,  is  the  derivation 
of  coefficients  for  iterative  methods  for  large  matrix  problems  Ax  =  b.  Soon  after  I  be¬ 
came  interested  in  this  kind  of  application,  I  realized  that  the  standard  way  of  treating 
matrices  in  these  problems — namely,  assuming  that  the  spectrum  defines  a  meaningful 
domain  for  conformal  mapping  or  other  approximation  methods — is  inappropriate  in 
some  cases  when  the  matrices  are  far  from  normal.  And  so  some  of  the  later  research 
covered  by  this  grant  has  been  in  a  new  area:  the  complex  analysis  and  numerical 
analysis  of  highly  non-normal  matrices.  This  work  is  summarized  in  Section  6. 

Most  of  the  projects  supported  by  this  grant  involved  extensive  computations  on  a 
network  of  Sun  workstations,  with  mouse  input  and  graphical  output.  I  believe  my 
early  computations  of  this  kind,  in  1985,  may  have  been  the  first  fully  interactive 
computations  of  conformal  maps  ever  carried  out  anywhere.  Five  year5  later,  so  much 
has  changed  technologically  that  that  statement  seems  astonishing. 
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2.  History  of  grant  AFOSR-87-0102 


3/1/85 

8/26/86 

12/1/86 

5/28/87 

7/22/87 

7/??/87 

12/1/87 

3/9/88 

6/28/88 

9/2/88 

12/1/88 

1/18/89 

2/26/90 


Proposal  submitted 
Announcement  of  initial  grant 
1st  year  starting  date  ($65,117) 

1st  year  Research  Progress  and  Forecast  Report  submitted 

Announcement  of  first  extension 

1st  year  Research  Summary  submitted 

2nd  year  starting  date  ($62,483) 

1st  year  Annual  Technical  Report  submitted 

2nd  year  Research  Progress  and  Forecast  Report  submitted 

Announcement  of  second  extension 

3rd  year  starting  date  ($66,400) 

2nd  year  Annual  Technical  Report  submitted 
Final  Technical  Report  submitted  (this  document) 


3.  Personnel 

As  the  Principal  Investigator,  I  was  supported  by  this  grant  for  three  months  per  year, 
enabling  me  to  reduce  my  teaching  load  by  one-third.  This  had  a  considerable  impact 
on  my  amount  of  time  available  for  research  on  numerical  conformal  mapping.  The 
work  was  distributed  throughout  the  year,  representing  somewhere  between  one- third 
and  one  half  of  my  research  output. 

A  graduate  student  named  Louis  Howell  worked  with  me  on  numerical  conformal  map¬ 
ping  throughout  the  three  years,  supported  primarily  by  this  grant.  Howell  finished 
his  PhD  thesis  last  month,  with  the  title  Computation  of  Conformal  Maps  by  Modified 
Schwarz-Christoffel  Transformations,  and  has  now  taken  a  job  at  Lawrence  Livermore 
National  Laboratory. 

Some  further  work  on  this  project  has  been  performed  with  the  collaboration  of  Noel 
Nachtigal,  a  third-year  graduate  student  in  our  group,  who  received  one  month  of 
support  from  this  grant,  and  Satish  Reddy,  another  third-year  graduate  student,  who 
did  not  receive  support  from  this  grant. 

In  the  final  year  this  grant  supported  a  two-month  visit  to  MIT  by  Lothar  Reichel 
(IBM  Bergen  Scientific  Centre  and  University  of  Kentucky),  which  led  to  several  of  the 
results  and  publications  described  below. 
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4.  Reports  and  publications 

The  following  reports  and  publications  were  supported  in  whole  or  in  part  by  this  grant. 

Copies  of  these  reports  will  be  sent  to  the  AFOSR  under  separate  cover.  Further  papers 

will  probably  be  written  later  by  Howell  concerning  the  topics  of  his  thesis,  [7]. 

[1]  F.  Dias,  A.  R.  Elcrat,  and  L.  N.  Trefethen,  “Ideal  jet  flow  in  two  dimensions,” 
Journal  of  Fluid  Mechanics  185  (1987),  275-288. 

[2]  L.  H.  Howell  and  L.  N.  Trefethen,  “A  modified  Schwarz-Christoffel  transformation 
for  highly  elongated  regions,”  SIAM  Journal  on  Scientific  and  Statistical  Comput¬ 
ing,  to  appear. 

[3]  L.  N.  Trefethen,  “Approximation  theory  and  numerical  linear  algebra,”  in  J.C. 
Mason  and  M.G.  Cox,  eds.,  Algorithms  for  Approximation  II,  to  appear. 

[4]  L.  N.  Trefethen,  “Schwarz-Christoffel  Mapping  iri  the  1980s,”  Numerical  Analysis 
Report  89-1,  Dept,  of  Mathematics,  M.I.T.,  January  1989. 

[5]  L.  N.  Trefethen,  “SCPACK  User’s  Guide,”  Numerical  Analysis  Report  89-2,  Dept, 
of  Mathematics,  M.I.T.,  January  1989. 

[6]  S.  C.  Reddy  and  L.  N.  Trefethen,  “Lax-stability  of  fully  discrete  spectral  meth¬ 
ods  via  stability  regions  and  pseudo-eigenvalues,”  in  Proc.  International  Conf.  on 
Spectral  and  High  Order  Methods,  to  appear. 

[7]  L.  H.  Howell,  Computation  of  Conformal  Maps  by  Modified  Schwarz-Christoffel 
Transformations,  PhD  thesis  and  Numerical  Analysis  Report  90-1,  Dept,  of  Math¬ 
ematics,  M.I.T.,  January  1990. 

[8]  N.  M.  Nachtigal,  S.  C.  Reddy,  and  L.  N.  Trefethen,  “How  fast  are  nonsymmetric 
matrix  iterations?”,  to  be  submitted  to  SIAM  Journal  on  Scientific  and  Statistical 
Computing. 

[9]  N.  M.  Nachtigal,  L.  Reichel,  and  L.  N.  Trefethen,  “A  hybrid  GMRES  algorithm  for 
nonsymmetric  matrix  iterations”,  to  be  submitted  to  SIAM  Journal  on  Scientific 
and  Statistical  Computing. 

[10]  L.  Reichel  and  L.  N.  Trefethen,  “Eigenvalues  and  pseudo-eigenvalues  of  Toeplitz 
and  block-Toeplitz  matrices,”  in  preparation. 

[11]  L.  Reichel,  “Newton  interpolation  in  Leja  points,”  to  appear  in  BIT. 

[12]  L.  Reichel,  “Fast  QR  decomposition  of  Vandermonde- like  matrices  and  polynomial 
least  squares  approximation,”  to  appear  in  SIAM  J.  Matrix  Anal.  Applies. 
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5.  Technical  results:  numerical  conformal  mapping 

5.1.  An  algorithm  for  computing  ideal  2D  jet  flows 

The  first  topic  of  research  supported  by  this  grant  was  the  computation  of  ideal  two- 
dimensional  free- streamline  jet  flows  by  conformal  mapping  (joint  work  with  Frederic 
Dias  and  Alan  Elcrat).  The  mathematical  formulation  goes  back  to  Helmholtz  in  1868, 
and  in  certain  situations,  these  idealized  models  can  yield  excellent  predictions  of  fluid 
behavior.  Potential  applications  include  design  of  airfoils,  hydrofoils,  and  propellers 
(analysis  of  cavitation),  and  various  other  hydrodynamical  problems  involving  weirs, 
blades  and  dams. 

Traditionally,  problems  of  this  kind  have  been  solved  by  use  of  the  hodograph  domain, 
but  that  method  requires  a  separate  analysis  of  every  flow  geometry,  the  calculation 
of  double  integrals  for  every  point  of  the  flow,  and  the  use  of  Riemann  surfaces  that 
may  be  intractable.  A  principal  contribution  of  our  work  was  the  observation  that  the 
introduction  of  a  modified  Schwarz-Christoffel  integral  makes  it  possible  to  dispense 
with  the  hodograph  domain  and  obtain  a  faster  solution  involving  a  single  integral, 
with  no  need  for  separate  analysis  of  each  geometry  or  the  use  of  Riemann  surfaces. 
We  have  implemented  the  resulting  algorithm  in  an  interactive  Fortran  package  called 
JET,  which  computes  the  ideal  jet  flow  through  any  polygonal  nozzle  to  arbitrary 
precision. 

This  work  has  been  published  in  [1].  An  example  of  a  jet  flow  calculated  by  the 
new  method  is  reproduced  below.  The  nozzle  shown  is  actually  a  polygon,  with  the 
quarter-circular  lip  approximated  by  16  straight  segments. 
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5.2.  An  algorithm  for  mapping  highly  elongated  polygons 

A  long-standing  problem  in  numerical  conformal  mapping  has  been  the  treatment 
of  elongated  domains.  For  well-understood  reasons,  the  conformal  map  of,  say,  a 
disk  or  half-plane  onto  a  rectangle  of  aspect  ratio  25  is  generally  impossible  even  to 
represent  in  floating-point  arithmetic,  let  alone  to  compute  with.  The  problem  lies 
in  the  mapping  itself,  not  in  the  algorithm  intended  to  compute  it.  Many  conformal 
mapping  algorithms  run  into  trouble  caused  by  this  so-called  “crowding  phenomenon,” 
and  my  own  SCPACK  is  no  exception.  Unfortunately,  many  applications  in  fluid 
dynamics  and  electronics  require  the  use  of  elongated  geometries. 

During  1987,  Louis  Howell  and  I  implemented  an  idea  for  getting  around  this  prob¬ 
lem  for  the  important  special  case  of  highly  elongated  polygons:  a  modified  Schwarz- 
ChristofFel  formula  based  on  the  idea  of  mapping  from  an  infinite  strip  rather  than  a 
disk  or  half-plane.  It  works  very  well,  enabling  us  to  map  regions  with  aspect  ratios  in 
the  tens  of  thousands.  Our  paper  [2]  on  this  subject  will  appear  in  the  next  issue  of 
SIAM  Journal  on  Scientific  and  Statistical  Computing. 

The  figure  below  illustrates  an  elongated  23-sided  polygon  that  has  been  mapped  to  a 
rectangle  by  the  new  method,  to  essentially  full  machine  precision.  (The  aspect  ratio 
of  the  rectangle  is  156.6241139.)  The  curves  plotted  are  conformal  images  of  horizontal 
and  vertical  lines  in  the  rectangle.  This  region  is  far  beyond  the  capabilities  of  any 
other  conformal  mapping  methods  that  I  know  of. 
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5.3.  Quadrature  methods  in  numerical  conformal  mapping 

Successful  implementation  of  the  Schwarz-Christoffel  transformation  and  other  related 
methods  of  numerical  conformal  mapping  depends  on  the  rapid  and  reliable  evaluation 
of  integrals  with  endpoint  singularities.  Many  methods  for  the  evaluation  of  such 
integrals  have  been  proposed  over  the  years.  Among  those  with  a  reasonable  claim 
to  be  close  to  optimal,  the  leading  candidates  are  compound  Gauss-Jacobi  quadrature 
(used  in  SCPACK)  and  various  more  fully  adaptive  methods  coupled  with  singularity 
removal  via  subtraction  (Kantorovich  and  Krylov)  or  change  of  variables  (Dias). 

These  quadrature  methods  have  been  compared  in  the  past,  but  not  from  a  modern 
point  of  view.  For  example,  a  PhD  thesis  by  E.-S.  Meyer  in  1979  devoted  many 
pages  to  a  comparison  of  methods  for  Schwarz-Christoffel  integrals,  but  ended  up 
rejecting  Gauss-Jacobi  methods  on  the  grounds  that  the  required  coefficients  are  hard 
to  obtain — which  is  simply  not  true  anymore,  thanks  to  software  such  as  GAUSSQ  by 
Golub  and  Welsch.  In  the  past  two  years,  therefore,  Howell  and  I  decided  to  carry 
out  an  extensive  and  more  up-to-date  comparison  of  various  integration  schemes  for 
Schwarz-Christoffel  problems.  We  believe  the  result  is,  relatively  speaking,  definitive. 

Our  conclusion  is  that  for  Schwarz-Christoffel  problems,  compound  Gauss-Jacobi  meth¬ 
ods  are  faster  than  other  integration  methods  by  a  factor  of  two  or  more.  This  factor 
can  be  thought  of  as  the  cost  of  adaptation,  for  the  Gauss-Jacobi  approach  avoids 
adaptation  by  using  a  priori  knowledge  ot  the  location  of  singularities. 

For  example,  the  lengths  of  the  sides  in  the  10-sided  polygon  below  can  be  evalu¬ 
ated  to  10-digit  accuracy  by  compound  Gauss-Jacobi  quadrature  with  112  integrand 
evaluations.  The  best  result  we  have  found  for  a  competing  method  is  538  integrand 
evaluations,  achieved  by  a  complex  version  of  the  adaptive  code  QUANC8  coupled  with 
a  change  of  variables  to  eliminate  singularities.  Lower  accuracy  specifications  narrow 
the  gap  somewhat,  but  compound  Gauss-Jacobi  quadrature  remains  the  winner. 

These  results  are  described  in  [7]. 
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5.4.  Divergence  of  Schwarz-Christoffei  algorithms 

In  practice,  SCPACK  nearly  always  constructs  conformal  maps  successfully — that  is, 
solves  the  Schwarz-Christoffei  parameter  problem  correctly.  The  same  is  true  of  the 
quite  different  algorithm  for  solving  the  parameter  problem  used  by  R.  T.  Davis  and 
others.  Nevertheless,  theorems  guaranteering  convergence  of  these  algorithms  have 
never  been  proved,  and  in  the  course  of  his  work  supported  by  this  grant,  Louis  Howell 
discovered  that  there  is  a  good  reason:  both  algorithms  diverge  in  some  cases.  These 
results  are  described  in  Howell’s  thesis,  [7j,  and  will  be  written  up  for  publication  later. 

I  will  illustrate  here  the  potential  divergence  of  SCPACK  only,  since  that  is  the  more 
surprising  result.  To  solve  the  parameter  problem,  SCPACK  adjusts  parameters  itera¬ 
tively  until  the  lengths  of  the  sides  of  the  polygon  are  correct.  This  amounts  to  solving 
a  system  of  nonlinear  equations,  and  to  do  this,  a  standard  piece  of  software  is  used 
(NS01A),  w'hich  guarantees  descent  of  an  appropriate  objective  function  at  every  step. 

The  example  below,  due  to  Howell,  explains  why  this  process  cannot  always  converge. 
Suppose  the  initial  guess  is  the  polygon  on  the  left  and  the  target  is  the  polygon  on 
the  right.  To  get  from  one  polygon  to  the  other  without  changing  any  of  the  angles, 
SCPACK  will  have  to  shorten  the  long  slits  a  great  deal  so  that  they  can  pass  around 
each  other.  In  doing  so,  the  side  lengths  will  temporarily  have  to  move  far  from  their 
correct  values.  In  other  words  there  is  a  potential  barrier  between  the  initial  guess  and 
the  target,  which  no  descent  method  based  on  a  side-length  formulation  can  cross. 


SCPACK  does  indeed  fail,  we  found,  when  applied  to  examples  like  this.  On  the  basis 
of  this  and  other  examples  we  have  devised,  we  now  believe  that  none  of  the  existing 
algorithms  for  Schwarz-Christoffei  mapping  are  globally  convergent  for  all  polygons. 

So  far,  despite  a  number  of  attempts,  Howell  and  I  have  not  been  able  to  devise  an 
algorithm  for  Schwarz-Christoffei  mapping  that  is  guaranteed  to  converge  in  all  cases. 
However,  thanks  to  examples  like  this  one  and  an  associated  theory  developed  in  his 
thesis,  we  believe  we  are  close  to  this  goal. 
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5.5.  Algorithms  for  mapping  circular  polygons 

Since  Schwarz  in  ISG9,  it  has  been  known  that  the  Srhwarz-Christoffel  integral  can  be 
generalized  to  an  ordinary  differential  equation  for  the  conformal  mapping  of  circular 
polygons:  planar  regions  bounded  by  straight  line  segments  and/or  circular  arcs.  Such 
regions  come  up  surprisingly  often  in  applications,  and  the  set  of  circular  polygons  also 
has  the  aesthetic  advantage  that  it  is  closed  under  Mobius  transformations. 

The  numerical  realization  of  this  Schwarzian  o.d.e.  is  quite  another  matter.  At  least 
twenty  independent  implementations  of  the  Schwarz- Christoffel  integral  have  been  at¬ 
tempted  over  the  years,  but  only  once  to  my  knowledge  has  the  Schwarzian  o.d.e.  been 
implemented  numerically  in  the  past — by  Bjorstad  and  Grosse  (SIAM  J.  Sci.  St  at. 
Comp.,  1987).  Unfortunately,  like  most  of  the  early  implementations  of  the  Schwarz- 
Christoffel  integral  itself,  their  attempt  met  with  only  partial  success.  The  problem  of 
mapping  circular  polygons  proves  to  be  highly  ill-conditioned,  at  least  in  its  standard 
formulation,  with  the  effect  that  the  convergence  of  the  Bjorstad/Grosse  implementa¬ 
tion  is  somewhat  problematical. 

In  his  thesis  [7],  Louis  Howell  made  substantial  progress  on  this  problem,  without 
solving  it  completely.  The  details  are  many  and  cannot  quickly  be  summarized,  but  here 
is  a  kind  of  engineering  summary:  the  conformal  mapping  of  circular  arc  polygons  is 
now  two  or  three  times  as  practicable  as  it  used  to  be,  thanks  to  Howell’s  innovations  in 
several  parts  of  the  problem.  The  definitive  treatment  of  the  circular  polygon  problem 
remains  to  be  carried  out,  but  I  expect  Howell’s  thesis  will  prove  the  foundation  for 
that  definitive  treatment  when  it  finally  appears. 

The  picture  below  shows  the  mapping  of  a  particular  circular  polygon  to  the  unit  disk. 
The  interior  curves  are  conformal  images  of  concentric  circles  and  radius  vectors  in  the 
disk. 
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5.6.  A  formula  for  mapping  highly  elongated  circular  polygons 

Section  5.2  mentioned  the  success  we  have  had  with  an  algorithm  for  mapping  highly 
elongated  polygons,  with  the  motivation  of  combating  the  troublesome  phenomenon 
of  "crowding.”  Combining  this  idea  with  the  material  of  Section  5.5  suggests  the 
possibility  of  new  algorithms  for  mapping  highly  elongated  circular  polygons. 

In  his  „hesis  [7],  Howell  derives  a  Schwarz-Christoffel  kind  of  formula  for  this  purpose. 
To  the  best  of  our  knowledge,  this  has  not  been  done  before.  This  formula  will  permit 
the  mapping  of  domains  like  the  one  sketched  below,  which  could  certainly  not  be 
handled  by  the  standard  Schwarzian  o.d.e.  in  floating-point  arithmetic.  We  have  not 
yet  implemented  this  idea  numerically. 
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5.7.  “Schwarz-Christoffel  Mapping  in  the  1980’s” 

I  have  not  made  much  progress  towards  completing  the  book  I  envision  entitled  Schwarz- 
Christoffel  Maps.  In  January  of  1989,  however,  I  delivered  an  invited  talk  at  the  the 
Complex  Analysis  session  of  the  American  Mathematical  Society  meeting  in  Phoenix 
on  ’‘Schwarz-Christoffel  Mapping  in  the  1980’s,”  and  out  of  that  talk  developed  a  col¬ 
lection  of  annotated  transparencies,  [4],  which  comes  surprisingly  close  to  a  sketch  of 
my  eventual  book — of  course,  with  the  details  missing.  This  is  not  a  paper  for  publi¬ 
cation,  but  I  have  circulated  it  widely  as  a  technical  report,  and  I  think  it  is  the  most 
comprehensive  survey  in  existence  of  modern  applications  of  Schwarz-Christoffel  maps. 


5.8.  Second  edition  of  SCPACK  User’s  Guide 

The  first  SCPACK  User’s  Guide  was  released  as  an  ICASE  Internal  Report  in  1983. 
SCPACK  has  proved  surprisingly  durable  since  then,  but  its  User’s  Guide  came  to 
seem  increasingly  old-fashioned.  In  January  of  1989  I  updated  the  User’s  Guide  in  a 
second  edition,  [5],  circulated  as  an  MIT  Numerical  Analysis  Report,  and  this  is  now 
the  standard  reference  for  users  of  SCPACK. 
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6.  Technical  results:  applications  of  complex 
analysis  in  scientific  computing 


6.1.  Matrix  iterative  algorithms  based  on  complex  approximation 

Since  the  1950s,  connections  have  been  seen  between  complex  analysis  and  the  con¬ 
vergence  of  iterative  algorithms  for  solving  large  systems  of  equations  Ax  =  b.  Midway 
through  the  course  of  the  research  supported  by  this  grant,  I  became  interested  in 
these  problems  because  of  the  possibility  of  using  numerical  conformal  maps  to  calcu¬ 
late  coefficients  for  such  iterations.  For  example,  suppose  the  spectrum  of  a  matrix  A 
is  known  to  lie  in  the  region  ft  below.  Then  parameters  for  an  effective  matrix  iteration 
can  be  derived  via  polynomial  interpolation  of  the  function  z~l  in  Fejer  points  along 
the  boundary  of  ft  — the  images  of  roots  of  unity  under  a  conformal  map  f(z)  of  the 
exterior  of  the  unit  disk  onto  the  exterior  of  ft.  The  figure  shows  64  such  Fejer  points, 
computed  with  a  version  of  SCPACK  modified  to  map  exteriors  of  polygons. 


My  interest  in  examples  like  this  led  ultimately  to  research  supported  by  this  grant  in 
several  directions.  The  first  was  the  preparation  of  the  report  [3],  which  is  a  summary 
of  various  connections,  both  old  and  new,  between  matrix  iterations  and  real  and 
complex  approximation  theory.  The  second  has  been  the  report  [8],  which  to  the  best 
of  my  knowledge  is  the  first  systematic  attempt  to  assess  and  compare  the  convergence 
properties  of  the  leading  parameter- free  matrix  iterations  for  nonsymmetric  matrix 
problems.  The  iterations  investigated  in  [8]  are  known  as  CGNR,  GMRES,  and  CGS. 
A  point  we  emphasize  in  this  paper  is  that  the  convergence  of  CGNR  depends  on 
the  singular  values  of  A ,  while  the  convergence  of  GMRES  and  CGS  depends  on  the 
eigenvalues — assuming  A  happens  to  be  normal  (i.e.,  having  orthogonal  eigenvectors) 
or  close  to  normal. 
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6.2.  Unreliability  of  eigenvalue-based  algorithms 

The  most  interesting  outgrowth  of  the  work  just  described,  however,  was  the  discovery 
that  it  is  a  mistake  to  pay  too  much  attention  to  the  eigenvalues  of  a  matrix  if  it  is  far 
from  normal.  In  particular,  the  standard  practice  of  basing  algorithms  on  eigenvalues 
or  estimates  of  eigenvalues  is  a  dangerous  one.  When  a  matrix  is  far  from  normal,  its 
eigenvectors  may  form  so  ill-conditioned  a  basis  that  the  eigenvalues  have  negligible 
importance  for  practical  computations.  What  matter  more  are  what  I  call  the  pseudo¬ 
eigenvalues:  those  numbers  A  with  the  property  that  A  is  an  eigenvalue  of  A  +  E  for 
some  perturbation  matrix  E  with  ||i?||  <  e. 

For  example,  let  A  be  the  200  x  200  upper  triangular  Toeplitz  matrix  whose  first  row 
is  (I,  3/4,  3/8,  3/16,...  ).  The  spectrum  of  A  is  just  the  singleton  {1},  which  might 
suggest  a  matrix  iteration  based  on  polynomials  of  the  form  (1  —  z)n.  However,  for 
small  e,  the  e-pseudo-spectrum  is  approximately  the  disk  \z  —  3/2|  <  1,  which  suggests 
the  quite  different  iteration  based  on  (1  —  \z)n .  The  figure  below  shows  that  the  first 
iteration  leads  to  geometric  divergence,  while  the  second  leads  to  geometric  convergence 
down  to  the  level  of  machine  epsilon. 


This  example  is  artificial,  but  matrices  equally  far  from  normal  arise  in  many  applica¬ 
tions — for  example  in  Gauss-Seidel  and  SOR  iterations,  in  the  modeling  of  convection- 
diffusion  problems,  and  in  the  solution  of  differential  equations  by  spectral  methods.  In 
all  of  these  situations  and  in  others,  our  work  shows  that  methods  of  complex  analysis 
may  be  useful,  but  only  if  they  are  based  on  a  set  of  pseudo-eigenvalues  rather  than 
the  set  of  exact  eigenvalues.  These  observations  are  being  published  in  [3]  and  [8].  In 
addition,  I  plan  quite  a  bit  of  further  work  in  this  area. 
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6.3.  Spectra  and  pseudo-spectra  of  Toeplitz  matrices 

The  example  on  the  last  page  illustrated  that  the  spectrum  of  a  Toeplitz  matrix  may  be 
very  different  from  its  pseudo-spectrum,  and  that  the  latter  may  matter  more  in  prac¬ 
tice.*  Toeplitz  matrices  arise  in  many  applications  in  differential  equations,  integral 
equations,  spline  approximation,  signal  processing,  and  other  fields.  During  Lothar 
Reichel’s  two-month  visit  to  MIT  in  1989,  supported  by  this  grant,  he  and  I  found  that 
Toeplitz  matrices  provide  an  interesting  family  of  matrices  with  which  to  study  the 
behavior  of  pseudo-eigenvalues.  For  example,  consider  the  100  x  100  Toeplitz  matrix 
with  1  on  the  super-diagonal,  1  on  the  sub-sub-diagonal,  and  0  elsewhere.  The  solid 
dots  in  the  figure  below,  lying  along  three  straight  spikes,  represent  the  spectrum  of 
this  matrix.  The  circles,  however,  show  what  happens  to  the  eigenvalues  when  random 
perturbations  of  size  t  =  10-4  are  added  to  the  matrix  elements.  These  are  examples 
of  e-pseudo-eigenvalues  for  a  particular  value  of  e,  and  in  any  computation  involving 
a  tolerance  of  this  order,  such  pseudo-eigenvalues  will  probably  have  more  practical 
significance  than  the  exact  ones. 


Reichel  and  I  are  preparing  a  paper,  [10],  containing  a  number  of  theorems  and  ex¬ 
amples  on  this  subject  of  spectra  and  pseudo-spectra  of  Toeplitz  matrices.  The  solid 
curves  in  the  figure  above  show  the  estimate  for  the  pseudo-spectrum  derived  in  this 
paper,  and  the  dashed  curves  show  what  happens  in  the  limit  N  — ►  oo,  where  we  get 
the  spectrum  of  the  associated  Toeplitz  operator. 

*A  Toeplitz  matrix  is  one  that  has  constant  entries  along  diagonals. 
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6.4.  A  hybrid  GMRES  algorithm  for  nonsymmetric  matrix  iterations 

Another  outgrowth  of  our  new  approach  to  eigenvalues  has  been  a  new  hybrid  iterative 
algorithm  for  large  nonsymmetric  systems  of  equations  Ax  =  b.  This  is  joint  work  with 
Reichel  and  my  student  Noel  Nachtigal.  In  investigating  existing  hybrid  algorithms, 
which  are  generally  based  on  estimating  eigenvalues  by  the  so-called  Arnoldi  process,  we 
found  that  breakdown  and  divergence  were  possible  because  of  the  same  unreliability 
of  the  use  of  eigenvalues  discussed  in  Section  6.2.  A  different  approach  based  on 
GMRES  rather  than  Arnoldi,  avoiding  eigenvalues  entirely,  turns  out  to  be  better 
both  theoretically  and  in  practice.  We  are  convinced  that  this  is  the  most  reliable  and 
also  the  most  natural  hybrid  iterative  algorithm  yet  developed. 

For  illustration,  the  figure  below  shows  the  performance  of  our  algorithm  when  applied 
to  a  Toeplitz  matrix  studied  in  a  paper  by  Joe  Grcar,  with  subdiagonal  entry  —1  and 
elements  1  on  the  main  diagonal  and  the  first  three  superdiagonals.  The  dimension  is 
N  =  1000,  and  the  figure  shows  the  error  as  a  function  of  work  for  restarted  GMRES 
and  our  hybrid  GMRES.  The  hybrid  algorithm  is  the  clear  winner. 

This  work  will  be  published  in  [9]. 
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6.5.  A  Lax-stability  criterion  based  on  complex  stability  regions 

Another  area  where  complex  analysis  appears  in  numerical  analysis  is  in  the  idea  of 
stability  regions  for  determining  stability  of  discretizations  of  ordinary  or  partial  differ¬ 
ential  equations.  In  particular,  suppose  a  time-dependent  partial  differential  equation 
is  discretized  in  space  by  finite  differences,  finite  elements,  or  spectral  methods,  then 
discretized  in  time  by  some  o.d.e.  formula  such  as  an  Adams- Bashforth  or  Runge-Kutta 
formula.  This  kind  of  separation  of  space  and  time  discretizations  is  called  the  Method 
of  Lines.  Will  the  result  be  stable?  The  conventional  wisdom  is  that  to  test  for  sta¬ 
bility,  one  should  test  whether  the  eigenvalues  of  the  spatial  discretization  operator 
lie  in  the  stability  region  for  the  time-stepping  formula.  It  has  long  been  recognized, 
however,  that  this  condition  is  necessary  but  not  sufficient  for  stability. 

My  student  Satish  Reddy  and  I  have  found  that  the  proper  way  to  analyze  such  prob¬ 
lems  is  once  again  to  replace  eigenvalues  by  pseudo-eigenvalues.  In  [6]  and  in  another 
paper  to  follow,  we  prove  theorems  to  the  following  effect: 

Theorem.  A  method  of  lines  discretization  is  stable  if  and  only  if  all  the 
e-pseudo-eigenvalues  lie  within  a  distance  0(e)  +  0(Af)  of  the  stability  region 
as  e  — » 0  and  At  — >  0. 

In  other  words,  roughly  speaking,  the  old  necessary  conditions  become  necessary  and 
sufficient,  if  eigenvalues  are  replaced  by  pseudo-eigenvalues. 

These  distinctions  have  important  consequences  for  practical  computations.  For  ex¬ 
ample,  the  picture  below  shows  a  solution  of  u<  =  ux  on  [—1, 1]  by  a  Legendre  spectral 
collocation  method.  Because  the  time  step  is  too  large,  a  instability  has  appeared  at  the 
boundary.  According  to  eigenvalue  analysis,  one  would  have  expected  this  computation 
to  be  stable,  but  the  pseudo-eigenvalues  tell  otherwise. 
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6.6.  Stable  algorithms  for  polynomial  interpolation  and  least-squares 

Finally,  I  will  briefly  mention  the  papers  [11]  and  [12]  by  Lothar  Reichel,  which  were 
completed  during  his  visit  to  MIT  supported  by  this  grant.  Both  concern  new  and  more 
stable  algorithms,  based  on  methods  of  complex  analysis,  for  problems  of  polynomial 
approximation.  In  [11],  Reichel  shows  that  whereas  polynomial  interpolation  on  a  real 
or  complex  set  can  be  highly  unstable  if  the  usual  Newton  interpolation  formula  is 
used,  the  instability  largely  vanishes  if  the  interpolation  points  are  chosen  and  ordered 
according  to  what  is  known  as  a  Leja  sequence.  In  [12],  he  investigates  a  new  and  more 
stable  algorithm  for  determining  polynomials  orthogonal  with  respect  to  an  inner  prod¬ 
uct,  which  leads  to  a  new  method  for  QR  decomposition  of  Vandermonde- like  matrices 
and  thence  to  an  improved  algorithm  for  polynomial  least-squares  approximation. 
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